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We study the noise characteristics of stochastic oscillations in protein number dynamics of simple 
genetic oscillatory systems. Using the three-component negative feedback transcription regulatory 
system called the repressilator as a prototypical example, we quantify the degree of fluctuations in 
oscillation periods and amplitudes, as well as the noise propagation along the regulatory cascade in 
the stable oscillation regime via dynamic Monte Carlo simulations. For the single protein-species 
level, the fluctuation in the oscillation amplitudes is found to be larger than that of the oscillation 
periods, the distributions of which are reasonably described by the WeibuU distribution and the 
Gaussian tail, respectively. Correlations between successive periods and between successive ampli- 
tudes, respectively, are measured to assess the noise propagation properties, which are found to 
decay faster for the amplitude than for the period. The local fluctuation property is also studied. 
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I. INTRODUCTION 

Molecular events driving cellular processes, such as the 
binding/unbinding of transcription factors to DNA sites, 
the transcription of mRNAs, the translation of proteins 
from the transcribed mRNAs, and their degradations, 
are fundamentally stochastic processes governed by prob- 
abilistic contacts between these molecules and intracel- 
lular diffusion [Ij. Therefore understanding the dynamic 
fluctuations in the genetic regulatory system arising from 
such stochasticity is crucial as they may provide a poten- 
tial obstacle in maintaining the robustness of the cellu- 
lar function and the homeostasis of the cellular behavior 
directly or indirectly. Study on stochastic fluctuations 
of cellular dynamics has been accelerated by the recent 
advances in synthetic biology [5, ^ and in experimental 
techniques using the fluorescent proteins [4], which has 
allowed us to monitor the fluctuating molecular activi- 
ties at the single-cell level in a controlled and designed 
manner. Theoretical framework with which stochastic 
dynamics of molecule number in time around the steady 
state can be analyzed has been developed and applied for 
simple gene expression systems How the fluctuation 
(or "noise") in one molecular component propagates to 
the downstream components in a simple, linear genetic 
network has also been studied @ . 

Most studies so far, however, concern fluctuations 
around the steady state, such as the variance in protein 
number around its temporally-averaged constant mean 
value. Meanwhile, there are only a limited number of 
experimental and theoretical studies concerning out-of- 
steady-state activities, such as oscillatory gene expres- 
sion Q- Therefore, their characteristics of dynamic fluc- 
tuations are largely unknown. To address this, here we 
study the characteristics of fluctuations in the non-steady 
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states, that of oscillatory systems in particular, by using 
a model genetic oscillatory system. Undoubtedly, oscil- 
latory system is one of the most important systems that 
never reach steady states, which underlie a number of im- 
portant cellular and physiological processes ranging from 
bacterial chemotaxis to mammalian circadian rhythms 
(sl-fToj. A number of synthetic gene circuits with a few 
components and relatively simple regulatory mechanisms 
have been constructed to realize oscillatory dynamics in 
protein numbers in vivo pdl - fTsj . Among these models, 
we consider the so-called repressilator (Fig. 1) [ll|, a sim- 
ple negative feedback circuit with three genetic compo- 
nents, for its simplicity in the regulatory processes in- 
volved, thus making it appropriate for a starting point 
in a systematic understanding of the noise properties of 
more complicated oscillatory systems. 

The repressilator system has been studied theoreti- 
cally in both determinisitc [13, [ill and stochastic [H, [l3| 
frameworks. Most of these works were concerned with 
the condition and the stability of sustained oscillations 
in the circuit with respect to various system parameters, 
such as DNA binding mechanisms, reaction rates, and 
plasmid numbers. A general conclusion from these stud- 
ies was that the stability and the coherence of the oscilla- 
tion depends strongly on the system parameters and that 
optimal conditions exist under which the system exhibits 
the most stable oscillations. However, even under the 
optimal conditions, the oscillation is never perfect, and 
"noise" is inevitable. It is, thus, necessary to know the 
general characteristics of the noise in oscillatory dynam- 
ics in such a stable oscillation regime, which is currently 
lacking. Our primary aim in this paper is to fill this gap 
by performing a detailed analysis of noisy oscillations in 
this system. 

Fluctuations in stable oscillatory dynamics can be clas- 
sified into two major classes: global and local noise. The 
former refers to fluctuations in the global variables of os- 
cillatory signals, such as the oscillation periods and am- 
plitudes, resulting from a series of temporally organized 
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FIG. 1: (Color online) Schematic illustration of the stochastic 
model of the repressilator assembled into a plasmid containing 
three genes, lad, tetR, and A-cI, and its elementary reaction 
steps. See the text for the details. 



molecular events. The latter is mainly due to momen- 
tary stochastic stepwise changes in the molecule num- 
ber, which occur locally around the overall oscillatory 
signal. Previous studies on dynamic fluctuations have 
mostly been concerned with the latter. For the oscilla- 
tory activities, however, the characteristics of global noise 
or fluctuation is equally, if not more, important, as both 
the robust control of periods and achievement of suffi- 
ciently strong and regular strength (amplitudes) of the 
signal would be required for an appropriate response of 
the cell under fluctuating and noisy condition. Thus, our 
focus in this work is toward the characterization of the 
global noise and its propagation, with a brief discussion 
on the local noise. 

This paper is organized as follows: In Sec. II, we de- 
scribe the repressilator model and the methods of numer- 
ical calculations. In Sec. Ill, we present the main results 
and discuss their interpretations and potential biological 
meanings. In the final section, we will summarize and 
conclude the paper. 



a typical value in the physiological condition following 
Ref . : The binding rate of repressor proteins to each 
operator site is = 1 nM^^s^^; the unbinding rate of 
repressor proteins from the singly-occupied and doubly- 
occupied operator sites are Ui — 224 s~^ and U2 — 9 s~^, 
respectively; the protein translation rate is = 0.167 s~^ 
per mRNA molecule; and the half-life (inverse decay rate) 
of an mRNA molecule is = 120 s. The remaining reac- 
tion steps are mRNA transcription and protein degrada- 
tion. We choose the transcription rate in the free (with- 
out inhibitor binding) state to be Ai = a s~^, and that in 
the inhibited state to be much slower as Ao = 10~^ x Ai. 
Finally, the half-life of a protein molecule is set to be 
Tp = Tm/I3 (Fig. 1). In this model setting, the two pa- 
rameters a and /3 control the overall stability of the oscil- 
latory behaviors, the characteristics of which have been 
studied in terms of deterministic modeling [111, [3 EBl ■ 

B. Stochastic Simulations: Gillespie Algorithm 

To fully realize the stochastic trajectories of protein 
number dynamics, we perform dynamic Monte-Carlo 
simulations of the reaction processes. We use the ex- 
act stochastic method first developed by Gillespie in 
1970s ^T^. This algorithm proceeds as follows: Given 
a list of reactions that can occur at the moment {Ri} 
with their reaction rates {ri}, the reaction occurring next 
is chosen randomly in proportion to the reaction rates. 
Then, the chosen reaction, say Rk, occurs, resulting in 
changes in the molecule numbers involved in that reac- 
tion. At the same time, the time elapses by 5, which is 
determined by a random number sampled from the expo- 
nential waiting time distribution PwiS) = exp(— r^J), 
following from the assumption that the reaction process 
is a Poisson process with a rate rfe. We repeat these pro- 
cedures for a sufficiently long period of time to obtain a 
sample trajectory of the given reaction system. 



II. MODEL AND METHODS 
A. Repressilator Model 

The repressilator [ll| is a plasmid fed into an E. 
coli, consisting of three transcriptional repressor genes 
(Fig. 1). The first gene is lad from E. coli itself, en- 
coding the transcriptional repressor protein Lad. Lad 
inhibits the transcription of the second gene tetR from 
the tetracycline-resistance transposon TnlO. TetR, the 
protein product of tetR, subsequently inhibits the tran- 
scription of the third gene cl originating from A-phage. 
Finally, its protein product CI inhibits lad expression, 
completing the tri-component negative feedback loop. 

Elementary reaction steps involved in this process 
include the binding /unbinding of repressor proteins, 
mRNA transcription, translation into their protein prod- 
uct, and degradations of mRNA and protein molecules. 
The reaction rate for each of these steps is chosen as 



C. Stability of Stochastic Oscillations 

We perform the stochastic simulation for repressila- 
tor dynamics by varying the parameters a = a/(3 and 
/3, each from 2^^^^ to 2^". Depending on the parameter 
values, the degree of stable oscillations exhibited by the 
stochastic dynamics varies [13, II3l- To address that, we 
consider the auto-corrclation function 

'^P>(t)'^Pi(*+A) 

as a function of the time lag A, where Pi(t) is the molec- 
ular number of protein i as a function of time, {x{t)) 
denotes the time average of x{t) over t, and ax{t) is the 
standard deviation thereof. For oscillatory signals, we 
have an oscillatory behavior in r(A). For imperfect os- 
cillatory signals, r(A) decays in amplitude with A while 
oscillating. The more regular the sustained oscillation. 
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FIG. 2: (Color online) (a) Stability diagram of the stochastic repressilator system. Plotted are the color-coded first-peak values 
of the auto-correlation function r(A) of the protein number time series obtained at each grid point, (b) Typical time-courses 
of the protein Lad number with time for three representative parameter sets indicated in (a): (T) (solid), (2) (dashed), and 
(3) (dotted), (c) The auto-correlation function of the protein number time series shown in (b), corresponding to the stability 
measures 0.83(1) (solid), 0.56(2) (dashed), and 0.14(3) (dotted). Numbers in parentheses denote the uncertainty in the last 
digit. 



the slower is the decay in the amplitude of r(A) with 
time; thus, the first peak in r(A) (A > 0) is higher. 
Therefore, we use the value of r(A) at the first peak as 
a stability measure for stochastic oscillations (Fig. 2). 



III. RESULTS 

A. Stability Diagram of the Stochastic 
Repressilator 

The so-called stability diagram has been used to dis- 
play whether the system exhibits sustained oscillations 
or not in the system parameter space [111, [TH| • For the 



D. Finding Oscillation Periods and Amplitudes 



Due to constantly fluctuating protein numbers in time, 
it is not trivial to pinpoint exact timings of the start and 
the end of each oscillatory cycle under stochastic oscil- 
lations. In this work, we devise an ad-hoc procedure to 
determine the oscillation periods and amplitudes from 
a noisy oscillatory signal by assigning top and bottom 
points of the oscillation in the time series. To this end, 
we look at protein number dynamics at a coarse-grained 
level by introducing a time window, inside which local 
variations are integrated out. The width of the time 
window is chosen to be order of the longest single re- 
action time scale, 10^ s, which we find suitable for our 
parameter choice, by comparing the result with manual 
assignments. We identify the maximum and the mini- 
mum points within each (non-overlapping) time window 
and use their time series to locate the oscillation peaks 
and valleys. Oscillation periods are then the time inter- 
vals between two successive peak positions, and ampli- 
tudes are the difference in protein numbers at the valley 
and the following peak position (Fig. 4(a)). Distributions 
of oscillation periods and amplitudes determined in this 
way are calculated from time series with a typical length 
10^ s. 
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FIG. 3: (Color online) Noise-induced sustained oscillations. 

(a) Deterministic protein number dynamics with the param- 
eter set {a = 2^,13 — 2®) shows damped oscillations which 
decay to zero amplitude in the long-time limit according to 
the linear stability analysis, (b) Stochastic protein number 
dynamics with the same parameter set as in (a) shows the ex- 
istence of noise-induced sustained oscillations in that regime, 
(c) The auto-correlation function of the time series shown in 

(b) , leading to the stability measure 0.75. 
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FIG. 4: (Color online) (a) Portion of the sample trajectory of the Lad protein number with time for (a, /3) = (1/2,1/4). 

(b) Distributions of the oscillation periods Pt{T) and (c) that of the oscillation amplitudes Pa (A) in the protein number 
oscillations obtained from the full simulation with the same parameter set in (a). The horizontal bin sizes are (b) 100 s and 

(c) 100 molecules. The inset of (b) shows Pt{T) in semi-logarithmic scale, indicating an initial exponential increase of Pt{T). 
The solid curve in (c) is a fit to the Weibull distribution, Eq. (3). 



stochastic dynamics, one cannot set a strict borderline 
between the oscillating and the non-oscillating states due 
to the noisy dynamics. However, it is meaningful to 
ask if the overall stability characteristics of the oscilla- 
tory dynamics is affected by the stochasticity. In gen- 
eral, stochasticity elevates the degree of nonlinearity in 
the dynamics, thereby providing an additional potential 
source of instability leading to oscillatory behaviors [l^ . 
In other cases, noise can quench the system's dynamics, 
thereby destroying the oscillatory activity [20j . 

Figure 2(a) shows the stability diagram correspond- 
ing to stochastic oscillation of the repressilator, plotting 
the stability measure, the first peak value of r(A), in 
(5,/3) space (Figs. 2(b,c)). There is a clear dark region 
near the center of the stability diagram with decreasing 
stability away from the center. This central region cor- 
responds to the oscillation parameter region obtained in 
the deterministic modeling [HI, [l^. There is, however, 
an additional dark region near {a « 2^,/3 « 2^) with as 
high stability measure as the central one, suggesting the 
existence of sustained oscillations due to the fluctuation- 
driven instability (Fig. 3), an example of the so-called 
deviant effects |21]. In this additional region, which is on 
the border of the sustained oscillation regime in the de- 
terministic dynamics, large amplitudal fluctuations are 
expected. This additional oscillatory region, however, 
although theoretically intriguing, may not be physiologi- 
cally relevant because it corresponds to extremely unsta- 
ble proteins with very short half-lives. In typical cells, 
the protein half-life is known to be typically much longer 
than that of mRNA, thus excluding large values of /3. 



In the following, we fix the parameters to be a = 1/2 
and /3 = 1/4 (indicated by x in Fig. 2(a)), well inside the 
sustained oscillation region with the stability measure of 
0.63, to study the fluctuation properties of the stochastic 
oscillations in the system. These parameter values are 
chosen to be close to those used specifically in Ref. [ll| , 
as they reproduce well the in-vivo repressilator dynam- 
ics. We also verified that in this regime the stochastic 
simulation reproduces well the on-average characteristics 
of the deterministic dynamics such as the mean period 
and the mean amplitude, allowing us to focus exclusively 
on the additional role of stochasticity. 



B. Single Protein-Species Level Fluctuation: 
Distribution of Oscillation Periods and Amplitudes 

Due to the topological symmetry of the repressilator 
system and our choice of identical reaction parameters 
for all three components, the statistical properties of the 
protein number dynamics of the three genes are identical. 
Thus, we focus on the dynamics of the first gene product, 
Lad. A portion of the sample trajectory of the Lad time 
series is shown in Fig. 4(a). The distribution of oscillation 
periods Pt{T) and that of amplitudes Pa{A) obtained 
from the full simulation are shown in Figs. 4(b,c). The 
Lad protein number oscillates in time with a mean period 
/It = 6612 s and a standard deviation ctt — 747. Thus, 
the coefficient of variation (CV) is 

CV{T) = ^ = 0.11. (2) 
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FIG. 5: (Color online) Noise propagation, (a) Sample time series of the protein number dynamics illustrating the correlation 
measurements, (b-g) Correlation between the oscillation periods (b-d) and amplitudes (e-g) along the regulatory cascade. Plot 
of the consecutive Tci vs. Thaci (b) and Txeta vs. Thaci (c); that of two consecutive TLaci's (d). Same plots are drawn for the 
amplitude correlation in (e-g). The Pearson correlation coefficients measured from the plots are r = 0.76 (b), 0.60 (c), 0.31 (d), 
0.59 (e), 0.29 (f), and 0.15 (g), respectively, (h, i) Plots of the correlation coefficients between oscillation periods (triangles) 
and between amplitudes (circles) as functions of the regulatory step number on a (h) linear scale and on a (i) semi-logarithmic 
scale. 



Looking closer, Pt{T) decays like a Gaussian for large T 
while the distribution is skewed leftward and increases as 
cx e^/^° for small T (inset of Fig. 4(b)). 

For the amplitude, the fluctuation is found to be 
stronger, as can readily be seen from the sample time 
series in Fig. 4(a). The mean amplitude is = 3659 
molecules, and the standard deviation is a a = 1588, cor- 
responding to a CV of CV{A) = 0.43, almost fourfold 
larger than the periods. This higher degree of amplitude 
fluctuation suggests the necessity for regulatory mech- 
anisms controlling the oscillation amplitude to build a 
more robust cellular oscillator. The distribution of oscil- 
lation amplitudes Pa{A) is found to be reasonably de- 
scribed by the Weibull distribution 

PA(A) = ^(^y"e-(^A)\ (3) 

with parameters k w 2.33 and A « 4118 (Fig. 4(c)). The 
Weibull distribution appears in a wide range of problems, 
from the extreme statistics [l^ to the return-interval 
distributions in natural and social phenomena [23] and 
fractures in disordered media [24f|, which may provide 



clues for the mechanisms underlying the amplitude fluc- 
tuations. 



C. Noise Propagation 

Noise propagation refers to how the fluctuations of up- 
stream reaction processes affect the fluctuation of the 
downstream processes, such as the effect of an mRNA 
number fluctuation on the corresponding protein number 
fluctuation ^] or that of a transcription factor number 
fluctuation on its regulated protein number fluctuation 
0. Most studies on noise propagation so far have been 
concerned with a simple linear cascade [6] or a simple 
autoregulated feedback structure [25]. 

Here, we are interested in how the global noise, such as 
the oscillation period and amplitude fluctuations, prop- 
agates down the regulatory cascade. To this end, we 
plot in Fig. 5 successive oscillation periods and ampli- 
tudes for the first three steps. For immediately following 
periods, the periods are well correlated (Pearson correla- 
tion coefficient r = 0.76; Fig. 5(b)). The correlation de- 
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creases monotonically (Figs. 5(c,d)), but not purely expo- 
nentially (Figs. 5(li,i)), leading to a correlation between 
two consecutive periods of the same protein number as 
high as r — 0.31. This means that the period fluctu- 
ation of a protein component does propagate down the 
regulatory cascade, and that the effect is strong enough 
not to be fully attenuated, not even after a complete 
feedback cycle. Thus, the oscillation dynamics of the re- 
pressilator exhibits a correlated pattern period-wise, even 
though the underlying reaction processes are assumed to 
be purely Poissonian. 

Meanwhile, the oscillation amplitude correlations 
(Figs. 5(e-g)) are found to decay purely exponen- 
tially with the regulatory step number n as r(n) sa 
exp(— n/n*), with the characteristic step n* w 1.8, a 
much faster decay than that of the oscillation periods 
(Figs. 5(h,i)). This result implies that the oscillation am- 
plitudes not only fluctuate more strongly individually but 
also are less tightly controlled along the system. Ampli- 
tude control in a stochastic oscillation system would be 
a more serious issue for robust functioning in this re- 
spect [26i] . 



D. Local Fluctuations 

Around the overall oscillatory change, there are con- 
stant ups and downs in the protein number, resulting 
from the inherently stochastic reactions (inset of Fig. 6). 
This part of the fluctuation, which we call the local noise, 
can also affect the system's function as it introduces 
many local maxima and minima within a global oscil- 
lation period and provides a source of error in sensing 
the system's state. The degree of local noise strongly 
depends on the individual reaction rates [l^ [13] and, in 
general, anti-correlates with the stability measure. Here, 
we are again interested in the characteristics of local noise 
in the stable oscillation regime. 

The general characteristics of the noise in protein num- 
ber dynamics are encoded in its spectral density S{f), the 
squared magnitude of its (complex) Fourier component 
with frequency /. In particular, the high-frequency be- 
havior of S{f) tells us about the local noise properties. 
In Fig. 6, we show the spectral density of the Lad pro- 
tein number. There we can see a clear peak located at 
/i « 1.5 X 10"'* Hz corresponding to the primary os- 
cillation frequency, in good agreement with the mean 
oscillation period in Fig. 4(b). This prominent peak is 
accompanied by additional peaks at integer multiples of 
/i with decreasing height values. The frequency region 
/ > 10^'^ Hz is dominated by the local noise, in which 
S{f) decreases algebraically with / as 

5(/)^r^ (4) 

implying that the protein number locally performs an 
uncorrelated random walk around the overall oscillatory 
trajectory. This implies that in the stable oscillation 



10^ 




10"^ 10"" 10"^ 10"^ 10"^ 10° 10'' 

Frequency f (Hz) 



FIG. 6: (Color online) Spectral density of the Lad protein 
number dynamics. The peak at /i « 1.5 x 10~* Hz is in good 
correspondence with the mean oscillation period. The data 
for / > 10"'^ are binned logarithmically for visual clarity. 
The slope of the dotted straight line is —2, drawn for the eye. 
(Inset) Zoom-up of the protein number trajectory highlighting 
the local fluctuation. 

regime, the global oscillatory dynamics of the repressi- 
lator is primarily brought about by a system-level long- 
range coordination of the reaction processes, with only a 
weak perturbation due to stochasticity. 

IV. SUMMARY AND DISCUSSION 

In this paper, we have studied the noise characteris- 
tics of stochastic (noisy) oscillations in the repressilator 
system for stable oscillation regime. Distinguishing the 
global and the local noises, we first found that even in 
this relatively loosely regulated system, the oscillation 
period is well controlled, fluctuating moderately as ap- 
proximately a Gaussian with a CV value of 0.11, whereas 
the amplitudes fluctuate much strongly, with a Weibull 
distribution and a CV value of 0.43. Likewise, the period 
fluctuation of a protein species is attenuated more slowly 
down the regulatory cascade, leading to a correlated be- 
havior in adjacent oscillation periods, a short-term collec- 
tive memory, within the feedback cycle. The amplitude 
fluctuation, however, propagates much more weakly, with 
a characteristic step number of 1.8, less than the feed- 
back cycle length. Finally, the local fluctuations entail an 
inverse-square spectral density at high frequencies, mean- 
ing that the protein number trajectory performs locally 
an uncorrelated random walk around the global oscilla- 
tory dynamics. Furthermore, we have checked that the 
overall noise characteristics do not change qualitatively 
when we run the repressilator for different parameter sets 
with varying degree of stability, or we consider other sim- 
ple oscillator circuits based on negative feedback, such as 
the 3-component cyclic circuit with two activators and 
one repressor, supporting the generality of the results. 

The repressilator is built with one of the simplest nega- 
tive feedback structures capable of sustained oscillations. 
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Oscillatory circuits functioning in real biological systems 
exhibit much more complex structures with multiple reg- 
ulatory components, which might have evolved to ensure 
robust control and fine-tuning of the oscillatory activ- 
ities. In this respect, we hope that the current work 
can be extended to systematic structural and dynamic 
studies towards the design principles of biological oscil- 
latory networks including the stochastic effects [27, 28]. 
Indeed, our result that the oscillation amplitude fluctu- 
ates strongly even in the most stable oscillatory regime 
in the simple circuit may imply the intrinsic tunability of 
the simple oscillatory circuit [29| on the one hand, but at 
the same time, it demands an improved oscillatory net- 
work design for robust and resilient oscillatory function. 
Areas of future research to this end would include, e.g., 
investigating the role of various regulatory mechanisms 
besides simple transcription-factor binding for noise con- 



trol, especially for the oscillation amplitude [26*1 , and the 
effect of a molecular network on the stochastic dynam- 
ics [l^ [13, [13, El . Further insights into the stochastic 
dynamics of the biological oscillatory system could be 
gained by studies of more complex synthetic or natural 
circuits [H, [13] and of the theory of competition-based 
ecological systems 0, [s^j • 
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